arXiv:1505.01963vl [math.NA] 8 May 2015 


Wave-type threshold dynamics 
and the hyperbolic mean curvature flow 


Elliott Ginder 

Research Institute for Electronic Science 
Hokkaido University 


Karel Svadlenka 
Department of Mathematics 
Kyoto University 


Abstract 

We introduce a method for computing interfacial motions governed by curvature dependent ac¬ 
celeration. Our method is a thresholding algorithm of the BMO-type which, instead of utilizing 
a diffusion process, thresholds evolution by the wave equation to obtain the desired interfacial dy¬ 
namics. We also develop the numerical method and present results of its application, including an 
investigation of the volume preserving and multiphase motions. 


1 Introduction 


In this report, we propose a numerical scheme for the computation of the so-called hyperbolic mean 
curvature flow 

-gp(t,s) =-K(t,s)v(t,s), 7(0, s) = 7o(s), — (0,s) = v 0 (s)u 0 . (1) 

Here 7 : [0, T) x I —> M 2 is a family of smooth curves, n denotes the curvature and v is the unit outer 
normal vector. 

This and similar types of hyperbolic problems were derived as models of oscillatory interface motions. 
For example, m derives a model equation based on Hamilton’s principle, considering stationary points 
of a geometrical action with local energy density consisting of kinetic and internal energy 


e = 


1 

2 



2 



If the initial velocity is normal to the interface, the velocity remains normal during the motion, which 
leads to the equation 


d 2 7 
dt 2 


= env — Ve. 


Employing the variational structure, the authors show that the flow is locally well-posed and study 
conditions under which blow-up occurs. In the case of graphs, they show uniqueness of solution satisfying 
certain entropy condition. 

Similar equation 


d 2 7 
dt 2 


■= kv — Ve, 


which is found to be related to a hyperbolic Monge-Ampere equation, is studied in B EES] with an 
emphasis on the development of singularities. The derivation of the equation is not given but its relation 
to equations for the motion of relativistic strings in Minkowski space is mentioned. 

On the other hand, regarding damped oscillations, which appear on the solid-liquid interface of some 
crystals during melting or crystallization, BM derived the model equation 


P ~dt + P v = (V> + ’ t P") K ~ G 


1 




(also see the references therein). Here, v, ^ are the normal velocity and normal acceleration, and ip, (3, p, F 
denote the interfacial energy, kinetic coefficient, effective density and a driving force for crystallization, 
respectively. 

Equation ([l]) is also addressed in [B], where local existence and relations for the evolution of geometrical 
quantities are shown. It is a special case of a model equation for the motion of bubbles obtained in [7], 

diL dp 

p — = —pv — anv + J — pu(\7 ■ u — v ■ Du ■ v) — u — . 

C LL CLL 

Here, ^ is the material derivative of the velocity vector, p denotes mass density, p is a factor correspond¬ 
ing to pressure, a is surface tension and / denotes additional sources of momentum. 

The numerical solution of the above hyperbolic mean curvature flows have been addressed only 
scarcely. Except for straightforward front tracking schemes, [5J develops a crystalline algorithm for the 
motion of closed convex polygonal curves. These methods cannot directly manage singularities, such as 
topological changes or the presence of junctions in the multiphase case, yet such problems can be resolved 
by adopting the level-set approach. A level-set method for curvature dependent accelerations based on 
the results of Sethian and Osher is presented in [)7J. In this method, a nonlinear ill-posed problem has to 
be solved and it is not clear how the ideas can be extended to the multiphase setting. 

Therefore, in this paper we aim at constructing a numerical scheme based on threshold dynamics 
of the BMO type. The BMO algorithm was presented in 12) and is based on the fact that frequent 
thresholding of the solution to the heat equation approximates the evolution of level sets according to 
the standard (parabolic) mean curvature flow. It is not obvious that such a method would be feasible 
for curvature accelerated motions. Nevertheless, investigating the formulae for the solution to the wave 
equation has revealed that level sets of the solutions whose initial condition is the signed distance to the 
interface, evolve with normal acceleration equal to their mean curvature, when the thresholding interval 
approaches zero. 

In the following, we will present and formally justify the resulting algorithm (which we call HBMO, 
the hyperbolic BMO algorithm) together with its device for propagating velocities over the redistancing 
step without the need for explicit computation of velocities. We then comment on potential extensions 
of our scheme to more general flows, such as those involving phase volume constraints or the motion of 
junctions. We also include several numerical results and numerical confirmations regarding the proposed 
method. 


2 The HBMO algorithm 

The proposed HBMO algorithm for numerical approximation of the motion 0 in the case of a planar 
closed curve reads as follows. 

Given: initial curve 70 , its normal velocity Vq, a final time T and a time step r = T/N. 

1. Extend vq in a suitable way to a neighborhood of 70 . 

2. For t £ [0,t] solve the initial value problem 

u tt (t,x) = Au(t,x), u(0,x) = d 0 (x), u t (0,x) —-v 0 (x), (2) 

where d o(x) is the signed distance function to jo- 

3. Define 71 as the zero level set of u(x,r). 

4. For n = 1,2,..., TV — 1 repeat 

(a) For t £ [0, r] solve the initial value problem 

x) = 2Au(t, x), u(0,x) = 2d n (x) - d n - 1 (x), u t (0,x) = 0, (3) 

where d n (x) is the signed distance function to j n . 
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(b) Define 7„+i as the zero level set of u(x,t). 


In this section, we explain the ideas behind the derivation of the above HBMO algorithm and provide a 
formal justification for its convergence. We divide the explanation into two parts. We begin by addressing 
the construction of the first curve 71 (steps 1 - 3 in the algorithm), which will clarify the reason why a 
threshold-type scheme for the wave equations leads to accelerations proportional to curvature. Then we 
treat further HBMO steps (step 4 in the algorithm), focusing on obtaining the propagation of interfacial 
velocities throughout the evolution. 


2.1 The first step of the HBMO algorithm 

Given an initial closed smooth planar curve 70 and its smooth initial normal velocity vo, we construct an 
approximation 71 of the curve’s evolution by HMCF at time r. This is done as follows: 

1 . Extend Vq to the neighborhood of 70 in R 2 as 

vo (x,y) = v 0 (x,y), 

where (x , y) is the orthogonal projection of (x, y) on 70, i.e., the nearest point to (x, y) on the curve 
7 o- 

2 . Solve the initial value problem 

u tt =c 2 Au, u( 0 ,x)=d(x), u t ( 0 ,x) =-v 0 (x), ( 4 ) 

where d[x) is the signed distance function to 70, and where we take c = 1 here, (but we keep the 
general coefficient c for later convenience). We remark that the velocity vo is defined only in a 
neighborhood of 70 but, since 70 is smooth and r can be taken small enough, this fact will not 
hinder the subsequent step due to the finite speed of propagation. 

3 . Define 71 as the zero level set of u(x,t). 

We now examine the evolution of the interface, making use of the explicit representation formula. 
The solution to © is equal to u(t,x ) = v(ct,x ), where v is the solution of 

v tt = Au, v( 0 ,x) = d(x), v t ( 0 , x) = -v 0 (x)/c. 


Here and in the sequel, x = (xi,X2) and y = (2/1,2/2) are points in R 2 . Using the Poisson formula for v 
and a change of variables, we find that the formula for the original solution u reads 


u(t , x ) 


1 f d(y) + \ 7 d{y) ■ {y - x) - tv 0 {y) 
2 -Kct J B (x,ct) sj&t 2 -\y — x \ 2 


( 5 ) 


In order to analyze the resulting motion, we take a point on the interface and rotate and translate 
the coordinate system, so that the point becomes the origin and the outer normal to the interface points 
in the direction of the x’2-axis. We assume that the interface is smooth in a neighborhood of the origin. 
Then the value of the signed distance function at any point y G B(x,ct) can be approximated by the 
following Taylor expansion (see mi provided that x in ([ 5 ]) is close to the origin and t is sufficiently small 
(depending on the smoothness of the interface): 

d(yi,y 2 ) = y 2 + ^ y\ + ^K Xl yf - ^n 2 y!y 2 + ^2 e “(?/)?/“- 


Here k is the curvature of the interface at the chosen point (now the origin) and the e a ’s are smooth 
functions. We remark that the error functions, denoted by e a (with a multiindex), will vary from place 
to place and are always assumed to be smooth and bounded functions of their variables. 

We now proceed to calculate the contribution of each term in the signed distance expansion to the 
solution of wave equation u(t,x). 
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For the first term y 2 we have 


u 1 (t, x) = 


1 


y-2 + 


2/i - xi 
2/2 - *2 


2ttc< JB(x,ct) \J c 2 t 2 — |2/ — a: 

1 r 


dy 


B( 0,1) c 
£2 


27TCf , 

2 71 " Jb( o,i) \/l — M 2 
x 2 ,1 , ' 2 ’ r 
27r 
^2 


2(cfz 2 + aj 2 ) - x 2 ^ 2 




0 ^0 




dz 


. d6 dr 


Here we use the change of variables 2 = {y — x) /ct to transform the domain of integration to a fixed ball 
B( 0,1). Also we use the fact that the function z 2 /sjl — \z\ 2 is odd with respect to the x 2 -axis, hence its 
integral over B( 0,1) vanishes. 


By a similar calculation, the second term \tiy\ gives 


l (t,x) = 


1 K 


y\ + 


22/1 

0 


2/1 - X 1 
2/2 - X 2 


2irct 2 J B(XtCt) y'c 2 * 2 -!y-a:|2 


K 

47T , 

K , 


r B(0 


3 c 2 t 2 zf + 

,D 


dz 


= ^(A 2 + ^)- 


dy 


The third term \n Xl y\ contributes to the solution as follows: 


2 /? 


\t,x) = 


(f 


2/i - *1 
2/2 - X 2 


127ret J B (x,ct) sjc 2 t 2 - |y - x| 2 

K Xl r 9c 2 t 2 Xizf + x\ 

12ti Jb( 0 , 1 ) ^/l “ PI 2 

Zo„2,2 


dz 


6 


(3c 2 r + xf). 


dy 


The fourth term —^n 2 y(y 2 yields 


J (t,x) = - 


* + ( 2 T) 


47rct JB(x,ct) \Jc 2 t 2 —]y — x| 2 


dy 


K 

47r 


3c 2 t 2 x 2 z 2 + 6c 2 f 2 xiZiz 2 + xjx 2 


B(0,1) 




dz 


-^(c^+x 2 ). 
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The velocity term is approximated by its Taylor series as follows: 


’(t,x) = 


1 f tv 0 (y) 

2nct JB(x,ct) \/c 2 t 2 — |y — x \ 2 
t r vq (x + ctz) 

2^ Jb( o,i) \/l — |2p 

t 

2n 


dy 


dz 


«o(0) + V^o(0) • (X + ctz ) + X)|a|>2 e c*(t, z){t, x) 


-B (0,1) 




dz 


w o(°) i + afr u o( 0 )tei + ^2 e a (t,x)(t,x) a , 
M >3 


since |^(0) = 0. 


Finally, the error term in the signed distance expansion can be evaluated as 


u e (t,x ) = 


1 


e a (y)y a + v ( e a {y)y a ) • (y - x) 


2lret 


| ^ 4 /fl(x,ct) \Jc^t 2 - |j/ — a:| 2 


dy 


1 v /' 

27F | ^ 4 ^(°.l) 


e Q (a; + ctz)(a: + cte) a + V {e a (y)y a ) \ y=x-\-ctz ’ Ctz 


V^ 


ch 


= 51 e a(.t,x){t,x) c 

|a|=4 


The solution to the wave equation with initial condition d in the neighborhood of the origin can thus 
be written in the following way: 


u(t,x) = x 2 + — (c 2 f 2 + xf) + /ta ’ 1 ‘ Cl (3c 2 f 2 + xf) — K X2 ( c 2 t 2 + xf) — rt u (a:, f) + u e (x , f) 
2 6 2 

, 2 . 


= a- 2 + )7(c 2 f 2 + £ 2 ) + 1 (3c 2 f 2 + x 2 ) - 2i-^(c 2 f 2 + x 2 ) 

2 o 2 

-V 0 ( 0 )t- |yy( 0 )tei + ^ e a (t,x){t,x) a . 

i «|>3 


K X2 (-„ 2 j .2 


( 6 ) 


From this result we can make the following observation. If the distance travelled by the interface 
in the normal direction (i.e., the ^-direction) after time r is denoted by do, then this distance can be 
calculated from the relation u(r, 0, do) = 0. In particular, 

0 = u(t, 0 , d 0 ) = d 0 + ^kc 2 t 2 - F S 0 k 2 c 2 t 2 - v 0 (0 )t + ^ e a (r, d 0 )(r, d 0 ) a . 

M >3 


This relation implies that, in terms of the order in r, the second order approximation of do is co(0)t — 
|kc 2 t 2 , since the error term is of order 0(r 3 + dj)). Therefore, the error term can be estimated by 0(r 3 ) 
and solving the above equation for do, we obtain 


do = 


a>o(0)r — \k,c 2 t 2 + 0(t 3 ) 


1 — ^k 2 c 2 t 2 


= i> 0 ( 0 )t - -kc 2 t 2 + 0 (r 3 ). 


(7) 


This means that the interface moves with initial velocity Vo and with an acceleration equal to —c 2 k, the 
c 2 -multiple of the curvature. Thus taking c = 1, we obtain the desired approximation of the interface 
7i ~ 7( T )- 


2.2 Further steps of HBMO 

If we want to apply the above scheme to all time steps, the velocities along the interface would need to be 
computed. This would severely complicate the numerical solution, even if the velocity field was known. 
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The difficulties inherent to such an approach are made clear, for example, when considering evolutions 
that involve topological changes (e.g., interfaces that split apart, or contact each other). On the other 
hand, if we reset the initial velocity to zero, the interface will always accelerate from zero velocity and 
the accumulated velocities will not propagate. Hence, we need to provide a modification of the above 
scheme which inherits the velocity from the previous time step and accelerates the interface depending 
on its curvature. 

In order to design a scheme that avoids computation of interface velocities, it is necessary to account 
for the position of the interface at both the previous and the current time step. Thus, it is natural to 
consider a vanishing initial velocity in the wave equation and an initial condition of the form 

u(0,x) = ad n -i(x) + bd n {x). (8) 

Here d n denotes the signed distance function to the present interface 7 „, d n -1 is the signed distance 
function to the interface at previous time step 7 n -i, and both a and b denote real numbers (see figure [l] 
below). 



Figure 1: The ^-coordinate system. 


In order to find a precise expansion of the signed distance function d n _i to the previous interface, it 
is necessary to calculate the change in the direction of the outer normal to the interface. Let us denote 
the point on the previous interface 7 „_i by w n _i and the curvature at that point by n n -\. Analogous 
notation is also used for the next interface 7 „. We are now considering n = 1, but we still write the 
symbol n since we will be able to choose a and b to propagate the interface’s velocity over each time step. 
According to (|6|, the solution of the wave equation @ with initial condition d n _i in the ^-coordinate 
system of figure I] reads 


u n _i(t,^) 


6 + f (c 2 t 2 + &) + ^(3c 2 f 2 + £ 2 ) - (c 2 f 2 + £) - v 0 (0)t - v 0xi mi 

M>3 


We want to compute the unit outer normal v n to the level set of u n -\ at the point (0,d ra _i) in the 
^-coordinate system. Hence, 

_ 7M n -i(r,Q,il n -i) 

|Vu„_i(t, 0, <S„_i)|’ 

where 

Vu„_ i(t,£) = ^ n _i^i + ^y^ L (c 2 r 2 + ^)-K 2 _ 1 ^i^ 2 -i; Oxi (0)f,l 

+Vc ^ e a (t,m,a) a - 

M>3 


«n-1 , 


(cV+£) 
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This yields 

Vu„_i(t,0,<5„_i) = 


(Kn-l)xi 2 2 /n\ i K n-1 2 I 

- 2 - C T ~ ^0 X1 (0)r, 1-— C T 


+ ( ^ e a (r,(5 n _i)(r,5 n _i)“, ^ e a (r, <5„-i)(t, 5„_i) a ) , (9) 

| a |>2 | a |>2 


and 


/ . 1x9 ( ty ^ 1 ) T ; wher e ^ = -^ 0 :ci (0)t + V e a (r, c5„_i)(t, <5„_i) 

V 1 + (^n) 2 H >2 


( 10 ) 


From § we can calculate the deviation angle of the normal in one time step, 

„ , n 1 + X]|a|>2 e «( T ’^-l)(T, 2 \ 

cos6» = (0,1) • v n = i = 1 T 0( t ~), 

IT S|a|>2 ®c*(x, 3n — l)(x, dn— l) Q 

since <5 n _i = O(t). Hence, we conclude that 9 = O(r). 

The x-coordinate and ^-coordinate are rotated by the vector v n and thus, 


6 = 


£2 — 


1 


\/ ! + M) 2 

1 

y/l + K) 2 


(xi + ^x 2 ) 

(-^Xi + x 2 ) T 


We can now write the signed distance d n -\ in terms of x, as follows (for simplicity we write u, instead of 

O'- 

d n -l(x) = ?2 T T -(Kra-l)^!^! — "h e «(£)£ Q 

M= 4 

= TT^f (-"*1 + * 2 ) T <5„-i T 2 (^ 2 ) On T ra 2 ) 2 + ^=r(xi T vx 2 f 
+ ^2) 2 (^f=!r T<5„-i) T ^ e a (x,<5„_i)(x, 

|a|=4 

= X 2 + S n -1 + - K n _ iX 2 + -( K n _ i) Xl Xi - - K 2 _ 1 X 2 X 2 

—r/Xi — ^-x 2 + Tvzi(—+ 2i/xix 2 + ^ 2 x 2 ) T ^^^(3i/x 2 x 2 + 3^ 2 XiX 2 ) — ^ 2^3 

-^=i(-z/xf - 2v 2 x\x 2 + 2vx\x\ + ^ 2 x 2 ) T v 2 x\x 2 - ^j^rjixi + vx 2 ) 2 5 n -i 

T ^2 e a (zA-i)(Mn-i)“ + 'Y^e i (x)v l . 

|a|=4 i> 3 


Because of (10), the last error term can be written as 


^ ' e a (r, 5 n _i, x)(t, 5 n _i) . 
M>3 


Computing the solution of the wave equation with initial condition d„_i and zero initial velocity, we 
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obtain 


K n-1 / 2.2 i 2\ i ( K rt—l)xi a: 'l / 0 2.2 


Un-l{t,x) — X 2 + S n -i-\ -^—(c t + %) + 


6 2 

„2\ i o..™ ™ i , ,2 { J2X2 i „2\ 


■(c¥ + ^) 


—i^xi — ^^2 + ^¥ 2 (ct + x 2 ) + 2 i/xix 2 + v 2 (cr1r + x 2 )\ 

+ — 6 l),ri [3l/X2(c 2 t 2 + X 2 ) + ?>V 2 Xi(c 2 t 2 + X%)\ — V 2 Xi(?>C 2 t 2 + X?) 


6 

: ^= i [— vxi(3c 2 t 2 + x'f) — |i/ 2 x 2 (c 2 f 2 + xf) + 2vx\{c i i 2 + x 2 ) + v 2 x 2 (3 c 2 t 2 + x 2 )] 


2,2 1 2 \ 

C t + XI' 

„2.2 , 2 


~ 2 f 2 +x 2 ) 
,2/„2.2 , 2\ 


~ 2(l+J 2 ) ^~l[ C ^ + X 1 + 2i/XiX 2 + +Xj)] 

+ e a (t,x){t,x) a + 5Z e ct( T j Jn-I, t, x)(t, 5 n -i) C 


|a|=4 


|a|>3 


where we have used the calculations of the integrals from Poisson formula in Section 2.1 The solution 
for for the initial condition ad n -1 + bd n evaluated at the point (t, x) = (r, 0, S n ) thus reads 


/ n r \ I aK ri—l + bn n 2 2 aK n—l S~ 

u n \T■, 0, o n ) = (a + b)0 n + ao„_i H---cr-- 


nr 2 2 
5 n C T 


I v~ r 1 K „-1 . 2 r2 . (Kn-lja,! r 2 2 
— ~2~O n — 9 — 0 n + - 2 - vS n C T 


'{ 

-^-[-^Sn^T 2 + V 2 8l] - 2^pp2)5n-l[c 2 T 2 (l + V 2 ) + I/ 2 <£]} 

+ 5^ e a (r, 5 n )(r, S n ) a + 5^ e Q (r, <5 n _i, <5 „)(t, <5„_i) a . 

|a:|=4 | O' | >3 

The distance S n traveled by the interface in the normal direction then satisfies u n (t 7 0, S n ) = 0, which 
can be written as 


0 = -5 3 av 2K ^+8 2 av 2K "-' 




2 1 2 

+5 n (a + b - “^-i +fcK - c 2 r 2 _ |„2 + | i/c 2 T 2(( Kn _ 1 ) ;i;i _|_ 

+a5 n -1 + aK "-i+ bK " c 2 r 2 _ hizi a( j ji _ lC 2 r 2 _(_ error terms, 

where the error terms are of order 0(t 3 + 8 3 _ x + S 3 ). Moreover, the quantities 5 n -1 and v are of order 
0(t), at worst. Since the lowest order approximation of the above equation is 

8 n (a + b) + aS n -i + aK ™- 1 + bl '™ c 2 r 2 = q, 

we see that S n must be of order 0(<5„_i+r 2 ), which is at most 0(t). Hence the S 2 and <5 3 terms (including 
their coefficients) are of order 0(r 4 ) or higher and can be neglected. In this way we obtain the solution 
5 n to the equation as 


Sn = ~ 


aS n -i + aK "-i +bK " c 2 r 2 + 0 (t 3 + ^_ 1 ) 

a + b- aKl ~ 1 2 +bKl c 2 t 2 - p 2 


_ —a X _ 

— a+b^" 1 


(IKn-i+bKn 2 2 


2 (a+6) 


c 2 t 2 + 0 (t 3 ). 


It is natural to choose a, 6 so that — = 1 and a + b = 1, i.e., a = —1, b = 2. Moreover, setting 

c 2 = 2 yields 

<5 n = (5„_i - (2 k„ - k„_i)t 2 + 0(r 3 ), (11) 

which is the desired relation between S n and S n -i. In particular, if we assume that the curvature is 
changing smoothly, relation (111 can be rewritten as follows, 

S n = 8 n -!- K n T 2 - (K n - K n -!)T 2 + 0 (t 3 ) 

= Sn-! ~ K n T 2 + 0 (t 3 ). 


(12) 





























To see that this is a correct approximation, consider the one-dimensional point-mass motion with 
initial velocity v 0 and acceleration —n(t), thus solving x"{t) = —«(f), x(0) = 0, cc'(0) = v 0 . The solution 
is 


x(t) = Vot — / / k (u)duds. 

Jo Jo 


(13) 


Our algorithm (12), including the initial step from the previous subsection, thus gives the approximation 

fc-i 

2 .»u. i -v n -u- 2 -- u . -vx- . -v-. ^cm- 2 -!- 2 


<5i = vqt - E 0 t 2 + 0(r 3 ), S 2 = v 0 t - ^ k 0 t 2 - Kir 2 + 0(2r 3 ),..., 5 k = v 0 r - ^k 0 t 2 - t 2 ^ m + 0(fcr 3 ). 
The corresponding total distance is 

fc —1 fc—i 

" (14) 


c(fcr) = fcuoT — k-KoT 2 — -EE Kj + (fcr) 2 0(r), 

2 i=ii=i 


which converges as r —>• 0 to the function (13) by the trapezoidal rule. The above analysis also shows 
that the accumulation of errors through time does not spoil the approximation of the interface, at least 
until the development of singularities. 

Remark. Regarding the 3D problem, assuming that the signed distance can be expanded as 

d(y) = 2/3 + \ K iVi + \ K 2 yl + H.O.T ., 

then a similar approach using the solution to the three-dimensional wave equation suggests the validity 
of our method in higher dimensions: 


u(t,x) = 


1 


47rt 2 

1 


V%) -(y-x)) dS(y) 


[ (y3 + X «1 y\ + x« 22/2 + «i2/i(yi - ah) + ^ 22 / 2 ( 3/2 - X 2 ) + {V3 - X 3 )) dS(y) 

JdB(x,t) K 1 2 J 

— / (tz 3 + X 3 + —Ki(tZi + Xi) 2 + —K 2 (tZ 2 + X 2 ) 2 

t7r 7aB(o,i) v 2 2 

+KitZi(tZi + X\) + K 2 tz 2 (tz 2 + X 2 ) + tzoj dS(z) 

r / (^3 + Ei(3t 2 zf +x\) -(- E 2 (3t 2 2:| +X%j) dS(z) 

t7r JdB( 0,1) v 2 2 ) 

— f f (x 3 + -nixl + -k 2 x\ + Eit 2 sin 2 0c°s 2 <j> + E 2 t 2 sin 2 flsin 2 (j>) smOdOdcj) 

vk Jo Jo V 2 2 2 2 / 


47T 

1 

47T 


2 • 2tt(x 3 + ^Kix 2 + + ^(«i + k 2 )t 2 ^ ■ n 


= X 3 + Ei(l 2 +xf) + ^ 2 (t 2 + x 2 2 ). 

This yields the distance traveled by interface in time t as 

a = — ^(«i + «2 )t 2 = -Ht 2 , 

confirming that our scheme can also be used with the three-dimensional wave equation, and likely in any 
higher dimension. 


3 Numerical Tests and Properties 

Thresholding dynamical algorithms have the computational advantage that there is no need to calculate 
curvatures, and singularities are implicitly handled by the partial differential equation. Moreover, they 
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are extremely simple to implement (here, we need only solve the wave equation). Thus, one can construct 
the corresponding numerical method in a number of ways-we will make use of standard finite differences, 
as well as minimizing movements for use in investigating volume preserving motions. 

Our approximation method for computing the HMCF ([I]) requires one to solve a wave equation, for 
which there are a number of schemes. Since we are considering an interfacial motion embedded in the 
solution to a wave equation, care needs to be taken when detecting the interface and constructing the 
signed distance functions. In particular, the precise location of the interface is needed. Without this 
information, errors arising from its approximation will propagate via the evolution of the signed distance 
functions. 

In the same vein, since the governing partial differential equation is hyperbolic, care must be taken 
in the numerical implementation of the method when constructing initial conditions. In particular, 
phase locations are encoded by vectors at the nodes of the discretized domain, and so the finite element 
assumptions can introduce initial interfaces that are not smooth. One approach to alleviating issues 
caused by this is to use one step of the original BMO algorithm to smooth the initial interface. This is 
easily done, since it only requires one to solve the vector valued heat equation for a small time. 

We will first perform a convergence test for an idealized version of our algorithm. This test assumes 
that no errors are introduced under spacial discretization. Our second test introduces the spacial dis¬ 
cretization, but computes signed distance functions utilizing an idealized representation of the interface. 
The final test removes this idealization and examines our algorithm’s order of convergence using standard 
finite differences. 


3.1 Idealized numerical convergence 

This section investigates the convergence rate of our algorithm for a simple test problem, under a slight 
idealization. For a circle evolving by ([I]) with initial normal velocity vq (we assume no < 0 for simplicity), 
one can compute the evolution of the circle’s radius r(t) analytically by solving 

r"(t) = -r(0) = r 0 , r'( 0) = v 0 . (15) 

r(t) 

Multiplying the equation by r'(t) and integrating we obtain 

r'{t) = 21og r + Ci, 


where C\ = 21ogro + Vq is determined by the initial conditions. Since 


dr 


•v/—2 log r + Ci 



we have 

e Gl ^ 2 \/f ei 'f ^ ~ log =t + C 2 , where C 2 = r 0 e^erf^-^L^ . 

In particular, C 2 = 0 for zero intial velocity. 

The solution can thus be expressed as 


r(t) = exp 




which yields a more simple form for zero initial velocity: 

r ( i ) = roe -H" 1 (v / f^)] 2 . 

Since the function erf -1 (a;) tends to infinity when x —> 1—, the radius vanishes in a finite time, which 
can be computed from the relation 

^e~ Gl/2 {t e + C 2 ) = 1. 
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After substituting values of integration constants, this gives 


te = r 0 




We remark that the extinction time for zero initial velocity is t e = r 0 


(see figure [i]). 


(16) 



Figure 2: Graph of the solution for rg = 1, vq = 0. 


Let us now calculate the numerical approximation due to our scheme. That is, we consider two circles 
of radii r 0 > r\, centered at the origin, which describe the initial conditions. We then compute the radius 
of the circle given by the 0-level set of the solution (after time t) to the problem 

utt = c 2 Au, u(0, x) = 2d\{x) — do 0*0) «t(0,x)=0. 

Here, c 2 = 2 and di is the signed distance function to the circle of radius r. ( , which can be written as 
di(x) = Ti — \x\. Hence, denoting r = 2r\ — ro, the initial condition is given by u(0,x) = r — |x| and the 
solution of the above initial value problem (sufficiently far away from the origin) is 

„r*r) = J-f r ~ _A f r - 2 \ x + ctz \ + l -wr ^ 

U[t ’ X) 2c-KtJ BM Jc 2 t 2 ~\y-x\ 2 y 2nJ B(0A) x/T^F 

Since the initial condition is radially symmetric and there is no interaction with a boundary, we can 
restrict the values of x to points on the x-axis of the form x = (£,0), where £ > 0. Then the above 
simplifies to 


u(t,£, 0) = r - 


1 j-i ^ 2 -k £ 2 +3ct p£ cos 9 + 2c 2 t 2 p 2 n ppp n 

2tt Jo Jo — p 2 + 2 ctpt; cos 9 + c 2 t 2 p 2 ^ 


The exact form of the solution to the above equation is difficult to find, since it involves elliptic integrals. 
Therefore, let us instead assume that t = r is small and expand the integrand as a Taylor series in time. 
To this end, let us denote 


a(t) = £ 2 + 3ctp cos 9 + 2 c 2 t 2 p 2 , b(t) = \/^+ 1 2ctp^cos^T+dH 2 /^ 2 , G(t ) 

and calculate the first and second derivatives of G at t = 0. 


a(t) 

Kty 


G( 0) 
G'(t) 
G\ 0) 
G"(0) 


2c£ 3 p cos 9 + 3c 2 £ 2 tp 2 (l + cos 2 9) + 8 c 3 £t 2 p 3 cos 9 + 2c 4 t 3 p 4 

W) 


2cpcos9, 

3c 2 £ 2 p 2 (l + cos 2 9) ■ £ 3 — 2c£^pcos9 • §£2cp£cos0 3c 2 2 2 

-^--- = — p sin 9 
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We can thus approximate the solution at points (£, 0) as 


u(t,£, 0) = p 2 (G(0) + G'(0)T + ±G"(0)T 2 + Q(T 3 )^d0dp 

P (t; + 2crp cos 0 + p 2 sin 2 9 + 0(r 3 )') cLO dp 
a/ 1 -p 2V ' 


To determine the position £ of the interface after time r, we need to solve the equation 

r-£- ^ t 2 = 0(t 3 ). 

The corresponding quadratic equation £ 2 — (r + 0(r 3 ))£ + |c 2 r 2 = 0 has two roots but since for vanishing 
times we want to recover the value £ = r, we take the corresponding root 



^ 2 


- + Vr 2 — 2 c 2 t 2 ^ + 0(r 3 ) = - + r-^ + 0(r 3 ) 


= r — -t 2 + 0(t 3 ). 
r 


Combining this result with our considerations regarding the first HBMO step, we obtain the following 
HBMO algorithm for the special case of a circle. 

1. Prescribe initial radius ro, initial velocity vq and time step r = t e /N , where f e is the extinction 
time (16) and N € N. 

2. Set 


n=r 0 - 


2 r n 


3. Repeat the following for n = 2,..., N: 

r n = 2r n „i - r n _ 2 - 


2 r„_i - r n _ 2 


(17) 


We test the resulting scheme for ro = 1 and vq = 0. The results are shown in the table below, where 
the error is taken at the half time to extinction, i.e., 

error = \r N/2 - r{t e /2 )| = \r N/2 - exp(-[er/ _1 (i)] 2 )|. 


We observe exactly linear convergence with respect to the time step r as predicted by (14), which shows 
that the proposed HBMO algorithm gives correct (i.e., convergent) result when the initial curve is a 
circle. 


division number N 

error at time t e /2 

convergence order 

10 

0.073199 

- 

10-4 1 

0.019332 

0.960 

10-4 2 

0.004884 

0.992 

10-4 3 

0.001224 

0.998 

o 

1 —1 

0.000306 

1.000 

10-4 5 

0.000076 

1.000 

10 • 4 B 

0.000019 

1.000 

10-4 7 

0.000004 

1.000 


Table 1: Results of the HBMO algorithm with discretized time only. 
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3.2 Convergence analysis under spacial discretization 

Our next tests concern the algorithm’s convergence under finite differencing. A circle of radius 3/10 is 
centered at (a :,y) = (1/2,1/2) with zero initial velocity, and the initial condition for our algorithm is 
therefore the signed distance function (taking positive values inside the circle) to this interface. The 
calculations are performed in the unit square under zero Neumann boundary conditions and we use grid 
points on a square lattice with spacing Ax = Ay = 1/(N — 1). A Delaunay triangulation of these points 
is performed, and we assume a linear interpolation of the grid values within each element. This allows us 
to define and determine the precise location of the interface. In particular, the location of the interface at 
each time step is given as a collection of line segments, each of which corresponds to the zero level set of 
the signed distance function across the elements. We remark that the geometry of the interface depends 
on the triangulation and that, as the interface shrinks, the relative resolution of the grid is a source of 
error. 

We take the HBMO time step as At = f e /2 9 (t e « 0.376) and the time step for the finite difference 
approximation of the solution to the wave equation is r = Af/(2 6 ). We then compute as follows: 

1. Construct the signed distance function do(x,y) to the initial circle. 

2. Set di(x, y) = do(x, y) (zero initial velocity), and repeat steps (3 — 5), for n = 1, ..., MAX. 

3. Compute the finite difference approximation to the following the wave equation: 

! Utt = Au in (0, At) x ft 

= 0 on (0, At) x dn (18) 

u(t = 0, x) = 2 d n — d n _ i, ut{t = 0, x) = 0 in fl. 


4. Compute r n+ 1 , the average distance from (1/2,1/2) of the line segments describing the zero level 
set of the solution to (18). 

5. Set d n +i(x, y) = y/ (1/2 - x) 2 + (1/2 - y) 2 - f n+ 1 . 

The error table is show in Table [ 2 J Here, L 2 -error designates to the square root of the value below, 

v e 

At Yl ~ 1 ) t ) - f nf , 

n—1 


where r(t) denotes the solution to (15), N e is the largest positive integer satisfying N e At < t e , and the 
value of f n is zero for any iteration after which the approximate solution’s radius becomes zero. 


Grid Resolution N 

L 2 -error 

convergence order 

16 

0.132320 

- 

32 

0.126222 

0.068067 

64 

0.113587 

0.152170 

128 

0.088544 

0.359340 

256 

0.051558 

0.780184 

512 

0.022104 

1.221893 

1024 

0.006538 

1.757309 


Table 2: Errors obtained for finite differences (using ideal distance functions). 


As our final numerical test, we replace the steps 4 and 5 of the previous algorithm with the construction 
of the exact signed distance function according to the exact configuration of the finite element space’s 
approximation of the interface. In both cases, we perform the numerical computations on various grid 
resolutions and investigate the corresponding convergence orders. 
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Grid Resolution N 

L 2 -error 

convergence order 

16 

0.131575 

- 

32 

0.124060 

0.084851 

64 

0.109666 

0.177921 

128 

0.084465 

0.376696 

256 

0.044312 

0.930634 

512 

0.021216 

1.062548 


Table 3: Errors obtained for finite differences. 


The error table is show in Table [3j where we use the same convention as in the previous test. We 
remark that explicit computation of the distance functions does not change the order of convergence. 

Figure [3] shows the evolution of a two phase HMCF, with zero initial velocity, as obtained by our 
algorithm. Contrast to parabolic mean curvature flow, we observe oscillations in the interfacial motion- 
the effect of acceleration is also evident in the motion of the phase region’s centre. 




Figure 3: Evolution by u' = —nn (time is from left to right). 


3.3 Multiphase motions 


Since the solution to the wave equation depends only on local information, one can formulate a multiphase 
algorithm which is akin to the original multiphase BMO algorithm in Q[2]. In particular, thresholding 
a collection of independently evolving wave equations, can formally be shown to compute multiphase 
interfacial motions with normal acceleration equal to curvature. Nevertheless, since we would also like to 
investigate multiphase volume preserving motions, we introduce a reformulation of this idea. In particular, 
by choosing a small time step At, we find a function u : fi — > R^ 1 solving the following vector valued 
wave equation: 


' u tt = A u 

ij = o 

rt*(0, x) = 0 

M t = 0,x) = 2zq - z e _ At 


in (0, At) x ft, 
on (0, At) x 9ft, 
in ft, 
in f2, 


(19) 


where N denotes the number of phases, ft is a smooth bounded domain in R rf , and the initial condition 
is defined by the following signed-distance interpolated vector field (see [4] and m) 


N i 

Z t( X ) = '^2iPiX{d\(x)>e/2} + - (2 + ^i( X )j PiX{-e/2<d t i (x)<e/2}- 


( 20 ) 


Here, \e is the characteristic function of the set E, the vector p i is the i th coordinate vector of a regular 
simplex in R w_1 for i = 1,..., N, and d\{x) denotes the signed distance function to the boundary of phase 
i at location x and time t: 


inf yedPf Ik ~ 2/11 
~ infap‘ Ik - 2/11 


if x elf, 
otherwise. 


( 21 ) 
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We remark that, when N = 2, equation (19) is scalar. At time At, in a process called thresholding , each 
phase region PA* is evolved as follows: 


Pf Kt = {x G SI : u(At,x) ■ p i > u(At,x) ■ p kl for all k € {1,..., A}}. 


( 22 ) 


The vector field Zq is then reconstructed using the boundaries of these sets and the initial condition 
for the wave equation is updated. The procedure is then repeated. The defining characteristic of the 


vector field (20) is that, upon choosing a location x positioned within a distance e/2 of an interface, say 
corresponding to the boundary of phase k. z e satisfies the property: 


sO) Pk = 40*0, 


(23) 


so that the profile in the direction normal to the interface is given by a signed distance function. This 


fact allows one to obtain the curvature of the interface by computing the Laplacian of (23). Moreover, 
as shown in m, this setting introduces a multiple well potential which allows one to express multiphase 
volume constrained motions. 

Our approximation method for computing multiphase interfacial dynamics thus repeats the PDE step 
(19) and the thresholding step (22). An example of a numerical result is shown in figure |4j Here, in 
contrast to parabolic mean curvature flow, we remark that interfaces are not smoothed, but oscillate. 






Figure 4: Multiphase evolution by HMCF (time is from left to right). 

As in the parabolic case [ITS], the method of minimizing movements can be used to inspect volume 
constrained interfacial motions as well. In particular, regarding volume-constrained motions, the min¬ 
imization formulation of our algorithm can allow one to formally include constraints, via penalization. 
Denoting the prescribed volume of region Pi by Ai , wave-type minimizing movements can be used under 
the following functional minimization (see, 0): 

— + ^ dx + - ^ |Aj - meas(P™)\ 2 . (24) 

' 6 i— 1 

Here and « n _2 are created from the wave equation’s initial conditions, e > 0 is a small penalty 

parameter and the areas corresponding to u are obtained from the sets 

P? = {% € H; u(x) ■ Pi > u(x ) • p J Vj}. 

Figure [5] shows a two phase volume-preserving motion. Contrast to the standard (parabolic) volume¬ 
preserving mean curvature flow, whose evolution approaches a circle, the interfacial dynamics considered 
here oscillate. Using the multiphase formulation of our algorithm, figure [6] below shows a six phase 
oscillatory motion. Compuations are performed within the unit square and we remark that volumes are 
preserved within an absolute error of 1(D 5 . 


P n(^) — / 
Jn 


\u 2lt n _i -)- Un— 

2h 2 


4 Conclusion 

We introduced a method for approximating interfacial motion governed by curvature dependent acceler¬ 
ation. Our method is a thresholding algorithm of the BMO-type which, instead of utilizing a diffusion 
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Figure 5: Volume preserving motion (time is from left to right). 



Figure 6: Evolution by multiphase volume preserving HMCF. 


process, thresholds an evolution by the wave equation. We obtained the desired interfacial dynamics by 
means of a combination of signed distance functions and the convergence order in time of our method was 
obtained. We also performed a numerical analysis of the algorithm and presented results of its computa¬ 
tional application, including investigations of multiphase and volume-preserving motions. The numerical 
results lead to the expectation of convergence for the volume-preserving threshold dynamical algorithm 
as well, which we will investigate in the future. 


5 Acknowledgements 

The work of the first author was supported by JSPS Grant Number 25800087 and the second by JSPS 
Grant Number 26870224. 


References 

[1] S. Essedoglu, S. Ruuth, R. Tsai, Diffusion generated motion using signed distance functions, Journal 
of Computational Physics. 229:4, pp. 1017 1042, 2010. 

[2] E. Ginder, K. Svadlenka. The discrete Morse flow for volume-controlled membrane motions. Advances 
in Mathematical Sciences and Applications, 22 No. 1 (2012), pp. 205-223. 

[3] E. Ginder, K. Svadlenka, A variational approach to a constrained hyperbolic free boundary problem, 
Nonlinear Analysis 71, pp. 1527 1537, 2009. 

[4] E. Ginder. A variational approach to volume-controlled evolutionary equations. Dissertation, 
Kanazawa University, 2012. 

[5] M. E. Gurtin, P. Podio-Guidugli, A hyperbolic theory for the evolution of plane curves, SIAM J. 
Math. Anal. 22, pp. 575-586, 1991. 

[6] C. He, D. Kong, K. Liu, Hyperbolic mean curvature flow, J. Differ. Equs. 246, pp. 373-390, 2009. 

[7] Myungjoo Kang, A level set approach for the motion of soap bubbles with curvature dependent velocity 
or acceleration. Dissertation, University of California Los Angeles, 1996. 


16 









[8] K. Kikuchi, S. Omata, A free boundary problem for a one dimensional hyperbolic equation, Advances 
in Mathematical Sciences and Applications 9, pp. 775-786, 1999. 

[9] D. Kong, K. Liu, Z. Wang, Hyperbolic mean curvature flow: Evolution of plane curves , Acta Math- 
ematica Scientia 29B(3), pp. 493-514, 2009. 

[10] D. Kong, Z. Wang, Formation of singularities in the motion of plane curves under hyperbolic mean 
curvature flow, J. Differential Equations 247, pp. 1694- 1719, 2009. 

[11] P. G. LeFLoch, K. Smoczyk, The hyperbolic mean curvature flow, J. Math. Pures Appl. 90, pp. 
591-614, 2008. 

[12] B. Merriman, J. K. Bence, S. J. Osher, Motion of multiple junctions: A level set approach, J. Comp. 
Phys. 112, pp. 334-363, 1994. 

[13] R.Z. Mohammad, K. Svadlenka, Multiphase volume-preserving interface motions via localized signed 
distance vector scheme, Discrete and Continuous Dynamical Systems - Series S, accepted, 2014. 

[14] H. G. Rotstein, S. Brandon, A. Novick-Cohen, Hyperbolic flow by mean curvature, Journal of Crystal 
Growth 198/199, pp. 1256-1261, 1999. 

[15] K. Svadlenka, E. Ginder, S. Omata, A variational method for multiphase volume-preserving interface 
motions, Journal of Computational and Applied Mathematics 257, pp. 157-179, 2014. 

[16] K. Svadlenka, S. Omata, Mathematical modelling of surface vibration with volume constraint and its 
analysis, Nonlinear Analysis 69, pp. 3202-3212, 2008. 

[17] A. Tachikawa, A variational approach to constructing weak solutions of semilinear hyperbolic systems, 
Advances in Mathematical Sciences and Applications 4, pp. 93-103, 1994. 


17 



